16 Spatially Continuous Data II

Introduction

NOTE: This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

In the previous practice you were introduced to the concept of fields/spatially continuous data.

For this practice you will need the following:

• This R markdown notebook.
• A file called “Walker_Lake.RData”

The file contains a dataframe with geocoded observations of a series of variables. The dataframe is in tidy format, that is, each observation is a row with its corresponding variables (coordinates, variables).

The Walker Lake dataset originally was used for teaching geostatistics in Isaaks and Srivastava’s An Introduction to Geostatistics.

The dataset contains a set of false coordinates X and Y, two quantitative variables V, U, and a factor variable T. The variables are generic, but you can think of them as measurements of pollutants.

In addition, the file includes two custom functions as follows:

1. point2voronoi(sp)

This is a function to obtain Voronoi polygons based on a set of points. It takes an argument sp (a SpatialPointsDataFrame) and calculates a set of Voronoi polygons. The value (output) of the function is a SpatialPolygonsDataFrame with the polygons.

2. kpointmeans(source_xy, z, target_xy, k, latlong)

This is a function to calculate k-point means. It takes a set of source coordinates (source_xy), that is, the coordinates of observations to be used for interpolation; a variable z to interpolate; a set of target coordinates (target_xy), the points to interpolate z; the number of nearest neighbors k; and a logical value to indicate whether the coordinates are latitude-longitude (the default is FALSE).

Learning objectives

In this practice, you will learn:

1. About intervals of confidence for predictions.
2. Using trend surface analysis as an interpolation tool.
3. The difference between accuracy and precision in interpolation.

Suggested reading

O’Sullivan D and Unwin D (2010) Geographic Information Analysis, 2nd Edition, Chapters 8 and 9. John Wiley & Sons: New Jersey.

Preliminaries

As usual, it is good practice to clear the working space to make sure that you do not have extraneous items there when you begin your work. The command in R to clear the workspace is rm (for “remove”), followed by a list of items to be removed. To clear the workspace from all objects, do the following:

rm(list = ls())

Note that ls() lists all objects currently on the worspace.

Load the libraries you will use in this activity:

library(tidyverse)

-- Attaching packages --------------------------------- tidyverse 1.2.1 --
v ggplot2 2.2.1     v purrr   0.2.4
v tibble  1.4.2     v dplyr   0.7.4
v tidyr   0.8.0     v stringr 1.2.0
v readr   1.1.1     v forcats 0.2.0
-- Conflicts ------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(spdep)

Loading required package: sp
Loading required package: Matrix

Attaching package: <U+393C><U+3E31>Matrix<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:tidyr<U+393C><U+3E32>:

    expand

Loading required package: spData

library(plotly)

Attaching package: <U+393C><U+3E31>plotly<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:ggplot2<U+393C><U+3E32>:

    last_plot

The following object is masked from <U+393C><U+3E31>package:stats<U+393C><U+3E32>:

    filter

The following object is masked from <U+393C><U+3E31>package:graphics<U+393C><U+3E32>:

    layout

library(deldir)

deldir 0.1-14

library(spatstat)

Loading required package: spatstat.data
Loading required package: nlme

Attaching package: <U+393C><U+3E31>nlme<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:dplyr<U+393C><U+3E32>:

    collapse

Loading required package: rpart

spatstat 1.55-0       (nickname: <U+393C><U+3E31>Stunned Mullet<U+393C><U+3E32>) 
For an introduction to spatstat, type <U+393C><U+3E31>beginner<U+393C><U+3E32> 

Begin by loading the data file:

load("Walker_Lake.RData")

You can verify the contents of the dataframe:

summary(Walker_Lake.t)

      ID                  X               Y               V         
 Length:470         Min.   :  8.0   Min.   :  8.0   Min.   :   0.0  
 Class :character   1st Qu.: 51.0   1st Qu.: 80.0   1st Qu.: 182.0  
 Mode  :character   Median : 89.0   Median :139.5   Median : 425.2  
                    Mean   :111.1   Mean   :141.3   Mean   : 435.4  
                    3rd Qu.:170.0   3rd Qu.:208.0   3rd Qu.: 644.4  
                    Max.   :251.0   Max.   :291.0   Max.   :1528.1  
                                                                    
       U           T      
 Min.   :   0.00   1: 45  
 1st Qu.:  83.95   2:425  
 Median : 335.00          
 Mean   : 613.27          
 3rd Qu.: 883.20          
 Max.   :5190.10          
 NA's   :195              

You can also check that the two functions have been loaded to your workspace.

Uncertainty in the predictions

In the preceding practice, you were introduced to three methods for interpolating spatially continuous data based on a sample of observations, namely Voronoi polygons, IDW, and k-point means.

These algorithms accomplish one of the objectives we typically have for spatially continuous data, since they provide point estimates. Recall that a field that was the outcome of a purely spatial process was defined as: \[ z_i = f(x_i, y_i) + \epsilon_i \]

A prediction of the field at a location where it was not measured was defined as a function of the estimated process and some random residual as follows: \[ \hat{z}_p = \hat{f}(x_p, y_p) + \hat{\epsilon}_p \] The first part of the prediction (\(\hat{f}(x_p, y_p)\)) is the point estimate of the prediction, whereas the second part (\(\hat{\epsilon}_p\)) is the random element.

The three methods seen previously do not, unfortunately, provide an estimate for the random element, so it is not possible to assess the uncertainty directly, which would depend on the distribution of the random term.

There are different ways in which some crude assessment of uncertainty could be attached to the point estimates. For example, a very simple approach could be to use the sample variance to calculate intervals of confidence. This could be done as follows.

We know that the sample variance describes the inherent variability in the distribution of values of a variable in a sample.

Consider for instance the distribution of the variable in the Walker Lake data:

ggplot(data = Walker_Lake.t, aes(V)) + geom_freqpoly()

Clearly, there are no values of the variable less than zero, and values in excess of 1,500 are rare.

The standard deviation of the sample is:

sd(Walker_Lake.t$V)

[1] 301.1554

The standard deviation is the average deviation from the mean. We could use this value to say that typical deviations from our point estimates are a function of this standard deviation (to what extent, it depends on the underlying distribution).

A problem with using this approach is that the distribution of the variable is not normal, and the distribution of \(\hat{\epsilon}_p\); the standard deviation is centered on the mean (meaning that it is a poor estimate for observations away from the mean); and in any case the standard deviation of the sample is too large for local point estimates if there is spatial pattern (since we know that the local mean will vary systematically).

There are other approaches to deal with non-normal variables, for instance Wilcox’s test, but some of the other limitations remain.

wilcox.test(Walker_Lake.t$V, conf.int = TRUE, conf.level = 0.95)

    Wilcoxon signed rank test with continuity correction

data:  Walker_Lake.t$V
V = 100580, p-value < 2.2e-16
alternative hypothesis: true location is not equal to 0
95 percent confidence interval:
 418.7 475.6
sample estimates:
(pseudo)median 
         447.2 

As an alternative, the local standard deviation could be used.

Consider the case of k-point means. The point estimate is based on the values of the k-nearest neighbors: \[ \hat{z}_p = \frac{\sum_{i=1}^n{w_{pi}z_i}}{\sum_{i=1}^n{w_{pi}}} \]

With: \[ w_{pi} = \bigg\{\begin{array}{ll} 1 & \text{if } i \text{ is one of } kth \text{ nearest neighbors of } p \text{ for a given }k \\ 0 & otherwise \\ \end{array} \]

The standard deviation could be calculated also based in the values of the k-nearest neighbors, meaning that it would be based on the local mean. Lets interpolate the field using the Walker Lake data. First create a target grid for interpolation, and extract the coordinates of observations:

target_xy = expand.grid(x = seq(0.5, 259.5, 2.2), y = seq(0.5, 299.5, 2.2))
source_xy = cbind(x = Walker_Lake.t$X, y = Walker_Lake.t$Y)

Interpolation using \(k=5\) neighbors:

kpoint.5 <- kpointmean(source_xy = source_xy, z = Walker_Lake.t$V, target_xy = target_xy, k = 5)

We can plot the interpolated field now. These are the interpolated values:

ggplot(data = kpoint.5, aes(x = x, y = y, fill = z)) +
  geom_tile() +
  scale_fill_distiller(palette = "OrRd", trans = "reverse") +
  coord_equal()

In addition, we can plot the local standard deviation:

ggplot(data = kpoint.5, aes(x = x, y = y, fill = sd)) +
  geom_tile() +
  scale_fill_distiller(palette = "OrRd", trans = "reverse") +
  coord_equal()

The local standard deviation indicates the typical deviation from the local mean. As expected, the standard deviation locally is usually lower than the standard deviation of the sample, and it tends to be larger for the tails, that is the locations where the values are rare.

The local standard deviation is a crude estimator of the uncertainty, again because we do not know the underlying distribution.

Other approaches based on bootstrapping could be implemented, but they are beyond the scope of the present discussion.

The issues assessing the level of uncertainty in the predictions with Voronoi polygons, IDW, and k-point means reflect the fact that these methods were not designed to deal explicitly with the random nature of predicting fields.

Other methods deal with this issue more naturally. We will revisit two estimation methods that were covered before, and see how they can be applied to spatial interpolation.

Trend surface analysis

Trend surface analysis is a form of multivariate regression that uses the coordinates of the observations to fit a surface to the data.

Lets see how this would work with a simulated example. We will begin by simulating a set of observations, beginning with the coordinates in the square unit region:

n <- 180
df <- data.frame(x = runif(n = n, min = 0, max = 1), 
                 y = runif(n = n, min = 0, max = 1))

Plot these locations:

ggplot(data = df, aes(x = x, y = y)) + 
  geom_point() + 
  coord_equal()

Lets now simulate a spatial process as follows:

df <- mutate(df, z = 0.5 + 0.3 * x + 0.7 * y + rnorm(n = n, mean = 0, sd = 0.1))

A 3D scatterplot can be useful to explore the data:

plot_ly(data = df, x = ~x, y = ~y, z = ~z, color = ~z) %>% 
  add_markers() %>% 
  layout(scene = list(
    aspectmode = "manual", aspectratio = list(x=1, y=1, z=1)))

We can fit a trend surface to the data as follows. In this case, the trend is linear:

trend.l <- lm(formula = z ~ x + y, data = df)
summary(trend.l)

Call:
lm(formula = z ~ x + y, data = df)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.194786 -0.064316  0.005245  0.058750  0.230605 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.50496    0.01975   25.57   <2e-16 ***
x            0.32199    0.02496   12.90   <2e-16 ***
y            0.69977    0.02461   28.44   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09631 on 177 degrees of freedom
Multiple R-squared:  0.8434,    Adjusted R-squared:  0.8417 
F-statistic: 476.7 on 2 and 177 DF,  p-value: < 2.2e-16

Given a trend surface model, we can estimate the value of the variable \(z\) at locations where it was not measured. Typically this is done by interpolating on a fine grid that can be used for vizualization or further analysis, as shown next.

We will begin by creating a grid for interpolation. We will call the coordinates x.p and y.p. These we generate by creating a sequence of values in the domain of the data, for instance in the [0,1] interval:

x.p <- seq(from = 0.0, to = 1.0, by = 0.05)
y.p <- seq(from = 0.0, to = 1.0, by = 0.05)

For prediction, we want all combinations of x.p and y.p, so we expand these two vectors into a grid, by means of the function expand.grid:

df.p <- expand.grid(x = x.p, y = y.p)

Notice that while x.p and y.p are of vectors of size 21, the dataframe df.p contains {r}21 * 21 observations.

Once that we have the coordinates for interpolation, the predict function can be used in conjunction with the results of the estimation. When invoking the function, we indicate that we wish to obtain as well the standard errors of the fitted values (se.fit = TRUE), as well as the interval of the predictions at a 95% level of confidence:

preds <- predict(trend.l, newdata = df.p, se.fit = TRUE, interval = "prediction", level = 0.95)

The interval of confidence of the predictions at the 95% level of confidence, in the form of the lower (lwr) and upper (upr) bounds:

summary(preds$fit)

      fit              lwr              upr       
 Min.   :0.5050   Min.   :0.3109   Min.   :0.699  
 1st Qu.:0.8325   1st Qu.:0.6404   1st Qu.:1.025  
 Median :1.0158   Median :0.8253   Median :1.206  
 Mean   :1.0158   Mean   :0.8241   Mean   :1.208  
 3rd Qu.:1.1992   3rd Qu.:1.0072   3rd Qu.:1.391  
 Max.   :1.5267   Max.   :1.3331   Max.   :1.720  

These values indicate that the predictions of \(z_p\) are, with 95% of confidence, in the following interval: \[ CI_{z_p} = [z_{p_{lwr}}, z_{p_{upr}}]. \]

A convenient way to visualize the results of the analysis above, that is, to inspect the trend surface and the interval of confidence of the predictions, is by means of a 3D plot as follows.

First create matrices with the point estimates of the trend surface (z.p), and the lower and upper bounds (z.p_l, z.p_u):

z.p <- matrix(data = preds$fit[,1], nrow = 21, ncol = 21, byrow = TRUE)
z.p_l <- matrix(data = preds$fit[,2], nrow = 21, ncol = 21, byrow = TRUE)
z.p_u <- matrix(data = preds$fit[,3], nrow = 21, ncol = 21, byrow = TRUE)

The plot is created using the coordinates used for interpolation (x.p and y.p) and the matrices with the point estimates z.p and the upper and lower bounds. The type of plot in the package plotly is a surface:

trend.plot <- plot_ly(x = ~x.p, y = ~y.p, z = ~z.p, 
        type = "surface", colors = "YlOrRd") %>% 
  add_surface(x = ~x.p, y = ~y.p, z = ~z.p_l, 
              opacity = 0.5, showscale = FALSE) %>%
  add_surface(x = ~x.p, y = ~y.p, z = ~z.p_u, 
              opacity = 0.5, showscale = FALSE) %>% 
  layout(scene = list(
    aspectmode = "manual", aspectratio = list(x = 1, y = 1, z = 1)))
trend.plot #%>% add_markers(data = df, x = ~x, y = ~y, z = ~z, color = ~residual)

In this way, we have not only an estimate of the underlying field, but also a measure of uncertainty for our predictions, since our estimated values are bound, with 95% confidence, between the lower and upper surfaces.

It is important to note that, although the confidence interval provides a measure of uncertainty, it does not provide an estimate of the prediction error \(\hat{\epsilon}_p\). This quantity cannot be calculated directly, because we do not know the true value of the field at location \(p\). We will revisit this point later.

For the time being, lets apply trend surface analysis to the Walker Lake data.

We will first calculate the polynomial terms of the coordinates, for instance to the 3rd degree (this can be done to any arbitrary degree, however keeping in mind the caveats discussed previously with respect to trend surface analysis):

Walker_Lake.t <- mutate(Walker_Lake.t,
                        X3 = X^3, X2Y = X^2 * Y, X2 = X^2, 
                        XY = X * Y,
                        Y2 = Y^2, XY2 = X * Y^2, Y3 = Y^3)

We can proceed to estimate the following models.

Linear trend surface model:

WL.trend1 <- lm(formula = V ~ X + Y, data = Walker_Lake.t)
summary(WL.trend1)

Call:
lm(formula = V ~ X + Y, data = Walker_Lake.t)

Residuals:
    Min      1Q  Median      3Q     Max 
-576.05 -241.79   -4.77  201.48 1055.98 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 589.5289    34.5547  17.061  < 2e-16 ***
X            -1.0082     0.1903  -5.297 1.82e-07 ***
Y            -0.2980     0.1727  -1.726   0.0851 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 292.1 on 467 degrees of freedom
Multiple R-squared:  0.06338,   Adjusted R-squared:  0.05937 
F-statistic:  15.8 on 2 and 467 DF,  p-value: 2.29e-07

Quadratic trend surface model:

WL.trend2 <- lm(formula = V ~ X2 + X + XY + Y + Y2, data = Walker_Lake.t)
summary(WL.trend2)

Call:
lm(formula = V ~ X2 + X + XY + Y + Y2, data = Walker_Lake.t)

Residuals:
    Min      1Q  Median      3Q     Max 
-579.37 -220.93    3.96  200.66  997.30 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 320.484250  70.904832   4.520 7.86e-06 ***
X2            0.001145   0.003191   0.359  0.71994    
X            -0.325993   0.910872  -0.358  0.72059    
XY           -0.006281   0.002244  -2.799  0.00534 ** 
Y             3.737955   0.722097   5.177 3.37e-07 ***
Y2           -0.011409   0.002188  -5.215 2.77e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 283.1 on 464 degrees of freedom
Multiple R-squared:  0.1258,    Adjusted R-squared:  0.1164 
F-statistic: 13.35 on 5 and 464 DF,  p-value: 3.571e-12

Cubic trend surface model:

WL.trend3 <- lm(formula = V ~ X3 + X2Y + X2 + X + XY + Y + Y2 + XY2 + Y3, 
                data = Walker_Lake.t)
summary(WL.trend3)

Call:
lm(formula = V ~ X3 + X2Y + X2 + X + XY + Y + Y2 + XY2 + Y3, 
    data = Walker_Lake.t)

Residuals:
    Min      1Q  Median      3Q     Max 
-564.19 -197.41    7.91  194.25  929.72 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -8.620e+00  1.227e+02  -0.070 0.944035    
X3           1.533e-04  4.806e-05   3.190 0.001522 ** 
X2Y          6.139e-05  3.909e-05   1.570 0.117000    
X2          -6.651e-02  1.838e-02  -3.618 0.000330 ***
X            9.172e+00  2.386e+00   3.844 0.000138 ***
XY          -4.420e-02  1.430e-02  -3.092 0.002110 ** 
Y            4.794e+00  2.040e+00   2.350 0.019220 *  
Y2          -1.806e-03  1.327e-02  -0.136 0.891822    
XY2          7.679e-05  2.956e-05   2.598 0.009669 ** 
Y3          -4.170e-05  2.819e-05  -1.479 0.139759    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 276.7 on 460 degrees of freedom
Multiple R-squared:  0.1719,    Adjusted R-squared:  0.1557 
F-statistic: 10.61 on 9 and 460 DF,  p-value: 5.381e-15

Inspection of the results of the three models above suggests that the cubic trend surface provides the best fit, with the highest adjusted coefficient of determination, even if the value is relatively low at approximately 0.16. Also, the cubic trend yields the smallest standard error, which implies that the intervals of confidence are tighter, and hence the degree of uncertainty is smaller.

Lets compare two of these models to see how well they fit the data.

First, we use create an interpolation grid. Summarize the information to ascertain the domain of the data:

summary(Walker_Lake.t[,2:3])

       X               Y        
 Min.   :  8.0   Min.   :  8.0  
 1st Qu.: 51.0   1st Qu.: 80.0  
 Median : 89.0   Median :139.5  
 Mean   :111.1   Mean   :141.3  
 3rd Qu.:170.0   3rd Qu.:208.0  
 Max.   :251.0   Max.   :291.0  

We can see that the spatial domain is in the range of [8,251] in X, and [8,291] in Y. Based on this information, we generate the following sequence that then is expanded into a grid for prediction:

X.p <- seq(from = 0.0, to = 255.0, by = 2.5)
Y.p <- seq(from = 0.0, to = 295.0, by = 2.5)
df.p <- expand.grid(X = X.p, Y = Y.p)

To this dataframe we add the polynomial terms:

df.p <- mutate(df.p, X3 = X^3, X2Y = X^2 * Y, X2 = X^2, 
               XY = X * Y, 
               Y2 = Y^2, XY2 = X * Y^2, Y3 = Y^3)

The interpolated quadratic surface is then obtained as:

WL.preds2 <- predict(WL.trend2, newdata = df.p, se.fit = TRUE, interval = "prediction", level = 0.95)

Whereas the interpolated cubic surface is obtained as:

WL.preds3 <- predict(WL.trend3, newdata = df.p, se.fit = TRUE, interval = "prediction", level = 0.95)

The predictions are transformed into matrices for ploting.

Quadratic trend surface and lower and upper bounds of the predictions:

z.p2 <- matrix(data = WL.preds2$fit[,1], nrow = 119, ncol = 103, byrow = TRUE)
z.p2_l <- matrix(data = WL.preds2$fit[,2], nrow = 119, ncol = 103, byrow = TRUE)
z.p2_u <- matrix(data = WL.preds2$fit[,3], nrow = 119, ncol = 103, byrow = TRUE)

Cubic trend surface and lower and upper bounds of the predictions:

z.p3 <- matrix(data = WL.preds3$fit[,1], nrow = 119, ncol = 103, byrow = TRUE)
z.p3_l <- matrix(data = WL.preds3$fit[,2], nrow = 119, ncol = 103, byrow = TRUE)
z.p3_u <- matrix(data = WL.preds3$fit[,3], nrow = 119, ncol = 103, byrow = TRUE)

This is the quadratic trend surface with its confidence interval of predictions:

WL.plot2 <- plot_ly(x = ~X.p, y = ~Y.p, z = ~z.p2, 
        type = "surface", colors = "YlOrRd") %>% 
  add_surface(x = ~X.p, y = ~Y.p, z = ~z.p2_l, 
              opacity = 0.5, showscale = FALSE) %>%
  add_surface(x = ~X.p, y = ~Y.p, z = ~z.p2_u, 
              opacity = 0.5, showscale = FALSE) %>% 
  layout(scene = list(
    aspectmode = "manual", aspectratio = list(x = 1, y = 1, z = 1)))
WL.plot2

And, this is the cubic trend surface with its confidence interval of predictions:

WL.plot3 <- plot_ly(x = ~X.p, y = ~Y.p, z = ~z.p3, 
        type = "surface", colors = "YlOrRd") %>% 
  add_surface(x = ~X.p, y = ~Y.p, z = ~z.p3_l, 
              opacity = 0.5, showscale = FALSE) %>%
  add_surface(x = ~X.p, y = ~Y.p, z = ~z.p3_u, 
              opacity = 0.5, showscale = FALSE) %>% 
  layout(scene = list(
    aspectmode = "manual", aspectratio = list(x = 1, y = 1, z = 1)))
WL.plot3

Alas, these models are not very reliable estimates of the underlying field. As can be seen from the plots, the confidence intervals are extremely wide, and in both cases include negative numbers in the lower bound. The uncertainty associated with these predictions is quite substantial.

Another question, however, is whether the point estimates are accurate. Lets add the observations to the plot:

WL.plot3 %>%
  add_markers(data = Walker_Lake.t, x = ~X, y = ~Y, z = ~V, 
              color = ~V, opacity = 0.7, showlegend = FALSE)

Alas, the trend surface does a mediocre job with the point estimates as well.

A possible reason for this is that the model failed to capture the spatial variability. To explore this, we will plot the residuals of the model, after labeling them as “positive” or “negative”:

Walker_Lake.t$residual3 <- ifelse(WL.trend3$residuals > 0, "Positive", "Negative")

ggplot(data = Walker_Lake.t, 
       aes(x = X, y = Y, color = residual3)) +
  geom_point() +
  coord_equal()

Visual inspection of the distribution of the residuals strongly suggests that they are not random. We can check this by means of Moran’s I coefficient, if we create a list of spatial weights as follows:

WL.listw <- nb2listw(knn2nb(knearneigh(as.matrix(Walker_Lake.t[,2:3]), k = 5)))

The results of the autocorrelation analysis of the residuals are:

moran.test(x = WL.trend3$residuals, listw = WL.listw)

    Moran I test under randomisation

data:  WL.trend3$residuals  
weights: WL.listw  

Moran I statistic standard deviate = 17.199, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.4633803457     -0.0021321962      0.0007325452 

Given the low p-value, we fail to reject the null hypothesis, and conclude, with a high level of confidence, that the residuals are not independent. This has important implications for spatial interpolation, as we will discuss in the following practice.

Accuracy and precision

Before concluding this practice, it is worthwhile to make the following distinction between accuracy and precision of the estimates.

Accuracy refers to how close the predicted values \(\hat{z}_p\) are to the true values of the field. Precision refers to how much uncertainty is associated with such predictions. Narrow intervals of confidence imply greater precision, whereas the opposite is true when the intervals of confidence are wide.

An example of these two properties is as shown in Figure 1.

Figure 1. Accuracy and precision

Panel a) in the figure represents a set of accurate points, since they are on average close to the mark. However, they are imprecise, given their variability. This is akin to a good point estimate that has wide confidence intervals.

Panel b) is a set of inaccurate and imprecise points.

Panel c) is a set of precise but inaccurate points.

Finally, Panel d) is a set of accurate and precise points.

Accuracy and precision are important criteria when assessing the quality of a predictive model.

This concludes Practice 15.